Theory of Boundary Effects in Invasion Percolation 
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We study the boundary effects in invasion percolation with and without trapping. We find that the 
presence of boundaries introduces a new set of surface critical exponents, as in the case of standard 
percolation. Numerical simulations show a fractal dimension, for the region of the percolating 
cluster near the boundary, remarkably different from the bulk one. In fact, on the surface we find 
a value of _D™'' = 1.65 ± 0.02 (for IP with trapping = 1.59 ± 0.03), compared with the bulk 



value of D" 



1.88 ± 0.02 {DI 



1.85 ± 0.02). We find a logarithmic cross-over from surface 



to bulk fractal properties, as one would expect from the finite-size theory of critical systems. The 
distribution of the quenched variables on the growing interface near the boundary self-organises 
into an asymptotic shape characterized by a discontinuity at a value Xc = 0.5, which coincides 
with the bulk critical threshold. The exponent t"'^^ of the boundary avalanche distribution for IP 
without trapping is r""*^ = 1.56 ± 0.05; this value is very near to the bulk one. Then we conclude 
that only the geometrical properties (fractal dimension) of the model are affected by the presence 
of a boundary, while other statistical and dynamical properties are unchanged. Furthermore, we 
are able to present a theoretical computation of the relevant critical exponents near the boundary. 
This analysis combines two recently introduced theoretical tools, the Fixed Scale Transformation 
(FST) and the Run Time Statistics (RTS), which are particularly suited for the study of irreversible 
self-organised growth models with quenched disorder. Our theoretical results are in rather good 
agreement with numerical data. 
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I. INTRODUCTION 

Recently, a large effort has been devoted to the study 
of Invasion Percolation (IP) [^]^. Compared to standard 
percolation [Q , IP has the advantage of describing the dy- 
namical evolution of the invading cluster as well as the 
final result. Furthermore, since a connectivity condition 
is naturally implemented in IP, its dynamics do not pro- 
duce extra, undesired, finite clusters, as happens instead 
in standard percolation Even if IP is more difficult 
to treat theoretically (because it presents a non-local, ex- 
tremal deterministic dynamics in a quenched disordered 
medium |§,||), it has been considered the paradigm of a 
large class of self-organised critical models . The Bak 
and Sneppen model for punctuated equilibrium and 
the Sneppen model for surface dynamics belong to 
this class. 

In the standard theory of critical phenomena, the role 
of boundaries has been intensively analysed |^ , and for 
many physical situations, ranging from Ising models to 
the more recent class of Self Organized Models [p|,p^, 
their presence produces a novel set of critical indices re- 
lated to the surface. The reason for the new behaviour 
consists in the lack of a microscopic layer in the system. 
This changes dramatically the microscopic interactions 



in the surface region of the system, yielding eventually 
to a macroscopically observable characteristic behaviour. 
The standard theory of finite size scaling of a thermody- 
namical system close to its critical point predicts in two 
dimensions ||ll| a logarithmic cross-over of the critical ex- 
ponents from the boundary to the bulk. Consequently, 
the effect of the boundary extends over the whole sys- 
tem. This is due to the strong correlations peculiar to a 
critical system. Such a study has already been done for 
standard percolation, and the results are available in the 
literature ||l^ , [l3|| , but no similar analysis has been per- 
formed for IP. Among the approaches applied to models 
with extremal dynamics, going from Mean Field treat- 
ment to a recently introduced technique called Run 
Time Statistics only the latter one, when com- 

bined with the Fixed Scale Transformation method , 
seems to be able to capture the subtle effects due to the 
presence of a boundary in the system. 

In this work we present numerical and theoretical ev- 
idence that a peculiar behaviour on the boundary takes 
place also in IP. Some of the results reported have al- 
ready been pubhshed ||l6|. Here, we would like to give 
a complete and detailed description of numerical results 
and of the derivation of the analytical results of the pre- 
vious Ref. |T6|. Moreover, we present new analytical and 
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numerical results, like the computation of the boundary 
avalanche exponent and the extension of our analysis to 
the case of IP with trapping, which has no analogue in the 
standard percolation model. In particular, the results for 
IP with trapping have no counterpart in standard perco- 
lation theory Q . 

From a qualitative point of view, the analogy between 
boundary effects in ordinary critical phenomena and IP 
can be easily understood by considering that boundary 
sites have fewer neighbors than bulk ones and hence fewer 
chances to invade a new region. Moreover, IP is a self- 
organised critical model and, as the evolution time tends 
to infinity, it can be considered in the same way as an or- 
dinary thermodynamical system when the temperature 
is tuned at the critical value Tc. The crossover between 
boundary and bulk fractal properties is shown by consid- 
ering intersections of the percolation cluster with straight 
lines parallel to the external boundary. This subset of the 
percolating cluster has a fractal dimension that varies 
with the distance from the boundary. Using some the- 
oretical tools introduced for the study of fractal growth 
processes, the Run Time Statistics (RTS) §,|4| and the 
Fixed Scale Transformation (FST) we are able to 
study analytically this behaviour, with an estimation of 
the boundary fractal dimension that is in rather good 
agreement with the numerical value. This is done for 
IP with and without trapping. In addition we study the 
avalanche dynamics near the boundary, for IP without 
trapping, and we compute both numerically and analyti- 
cally, by using the RTS and FST schemes, the boundary 
avalanche exponent r'*"''. 

Our results are presented in the following order. In 
section II, we present the definition of the model and a 
review of the numerical data. In section III, we describe 
the concepts underlying RTS and FST. In section IV we 
apply these methods to the computation of the boundary 
fractal dimension. In section V we compute the boundary 
avalanche exponent. In the last part we give a summary 
of the main topics. Appendix A is devoted to the deriva- 
tion of the RTS equations. 

II. THE INVASION PERCOLATION MODEL 

IP was introduced more than 10 years ago Q in order 
to describe the slow capillary displacement of a fluid (e.g. 
oil), the defender^ from a random porous medium due 
to another immiscible invading fluid (e.g. water), the 
invader. 

In general two cases are studied: 1) the medium is 
filled with an incompressible defender (e.g. oil), which is 
immiscible with the invader fluid; 2) the medium is filled 
with a defender with an infinite compressibility. In the 
former case the invader may trap regions of the defender: 
e.g. as the water advances, it can completely surround 
regions of the oil. These regions become disconnected 
from the other bonds occupied by the defender and, due 



to incompressibility, they become forbidden to the in- 
vader. This trapping effect lowers the fractal dimension 
of the percolating invader cluster. From an experimental 
point of view, trapping is connected to the phenomenon 
of "residual oil" , which is a great economic problem in 
the oil industry Q . 

The random medium is represented by a network of 
bonds corresponding to the throats connecting the pores 
of the medium. Let us assume, now, that the invader 
begins to displace the defender. Under the condition of 
a low and constant flow rate, the interface can be con- 
sidered to move one step at time, by invading the throat 
with the smallest section, i.e. the throat where there is 
the largest capillary force One can mimic this be- 
haviour by assigning a random section Xi (here we take 
an uniform distribution in [0,1]) to each bond i of the 
medium. The invading cluster evolves by occupying the 
bond with the smallest Xi on its perimeter. This is what 
is called a deterministic extremal dynamics. 

To study the behaviour at the boundary of this model, 
we performed some numerical simulations in the system 
shown in FIG.^ representing a sample of a two dimen- 
sional square lattice. To study the effect of only one 
boundary (e.g. the left one), we ensured isolation from 
the other one. To obtain this, we choose a lattice with 
size HL X 2L where H = 2,3,4, and the initial invader 
cluster is composed of the first L bonds of the bottom 
line, starting from the left boundary. The simulation 
stops when the cluster percolates the system, i.e. when 
the growth reaches the top of the sample. 

In FIG.|| a typical realization of this process is shown. 
The region of interest is the bottom-left one in 
where we can assume that the region is "frozen" with 
respect to the invasion process, i.e. the asymptotic fractal 
properties of the percolating cluster are well defined. For 
each value of the system size L we collected a set of 10'^ 
different realizations. In the region where statistics is 
collected, we study the fractal dimension of the sets of 
points obtained by intersecting the percolating cluster 
with lines parallel to the boundary, at a distance z from 
it. In this way, we are able to follow the cross-over of the 
fractal dimension of the cluster from the boundary to the 
bulk region. A standard box-counting procedure is used 
to compute the fractal dimension of the intersections. 
The behaviour of the fractal dimension d{z/L) of the 
intersections as a function of the normalized distance z/L 
from the left boundary is presented in FIG.|^ for L = 
256, H = 4 (i.e. 512 x 1024). In Table | we present 
the values of the boundary fractal dimension d"""^ for 
different system sizes L and different values of H. For 
the largest simulation L = 256, H ^ 4 (i.e. 512 x 1024), 
we obtain the result that the fractal dimension of this 
subset of the cluster passes from d""'' — 0.65±0.02 on the 
boundary to d''"' = 0.88 ± 0.02 in the bulk (at a distance 
z/L ~ 0.4 from the boundary), where d''"' represents the 
fractal dimension of the intersection far away from the 
boundary. A similar behaviour holds for smaller sizes L 
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and H as well. Since the dimension of the intersection 
set d obeys d = D — \ where D is the fractal dimension 
of the cluster, the last result is in agreement with the 
known value of D ~ 1.89. 

In order to explain such a slow crossover from d""*" 
to c?''"' we assumed that the number of occupied sites 
N{z, L) at a distance z from the boundary follows the 
finite-size scaling law: 

N{z,L)^L''''"f{z/L) (1) 

where one has f{z/L) - for z«L and 

f{z/L) — const for z >> L. Then in the first region we 
should have: 

d{z) = rf'^"'' + (d''"' - d''''')log{z)/log{L) (2) 

To test this scaling hypothesis we collapsed the curves 
relative to different L by plotting [d{z) — c?'^"''] log{L) as 
a function of z. The result depicted in FIG.^ shows a 
rather good collapse in the small z region. A similar be- 
haviour is found for IP with site trapping. In order to 
implement site trapping in our simulations, after each 
growth step a fictitious Laplacian field cf) is relaxed on 
the growing structure, with the following boundary con- 
ditions: = on the bottom boundary, the left bound- 
ary and the invading cluster, while ^ = 1 on the top 
boundary. In this way, all the bonds in a closed, trapped 
region are characterized by = 0. Then it is possible to 
recognize trapped bonds and to eliminate them from the 
list of bonds allowed to grow at the next step. Obviously, 
in this case the numerical simulations need much more 
time to be performed and we have been able to collect a 
smaller, but still significative, statistics with respect to 
IP without trapping (lO'^ clusters for L = 64, 2 x 10^ 
clusters for L = 128 and 10^ clusters for L = 256 each 
one for iJ = 2,3,4). In FIG. | we show the behaviour of 
the intersection dimension dtr{z/L) versus the normal- 
ized distance from the boundary, each simulation is for 
a values of H = 4. The fractal dimension c?t,. is com- 
puted on samples 512 x 1024 (i.e. L = 256 and H — 4) 
and passes from df"'' = 0.59 ± 0.03 on the boimdary to 
d'lf = D\f - 1 = 0.85 ± 0.03 in the bulk, which is agree- 
ment with the known value Df ^ 1.82 for site trapping 
101 . The data shown in Table || exhibit the same slow 
logarithmic cross-over found for IP without trapping. 

Other important quantities characterizing the dynam- 
ical properties of IP are the average distribution of 
quenched variables on the perimeter, called histogram 
^((a;), which gives an evidence of the self-organised na- 
ture of the model, and the avalanche-size distribution in 
the asymptotic critical state Q{s]Xc)^ where Xc is the 
"self-critical" threshold of the model. 

Let us start with the study of the histogram for IP 
without trapping. It is known |]l^ that for the bulk IP, 
the histogram distribution evolves in time from the initial 
flat shape, and self-organises into a step function with a 
discontinuity at a critical value a;^"' which depends on 



the details of the model and on the embedding dimen- 
sion ||l| and coincides with the critical threshold of the 
classical percolation in the same kind of lattice. For the 
two dimensional square lattice one has x^"' = \. Our 
simulations show that the distribution of the XiH on the 
boundary self-organises into a theta function and the 
critical threshold is again x^"'. A comparison between 
the bulk histogram and the boundary one is shown in 
FIG.^ It is not surprising to find a similar behaviour, 
because the value of the boundary critical threshold is 
dependent on the dynamical evolution of the whole per- 
colating cluster. Since for bulk IP the trapping does not 
affect the histogram distribution ||l[], the introduction of 
the trapping does not modify the above result. 

Another important quantity describing the dynam- 
ics of the model is the critical avalanche-size distribu- 
tion Q(s;Xc)- An avalanche is a sequence of elementary 
growths events causally and geometrically connected to 
a first one, which is called the initiator of the avalanche. 
That is, if one consider an event of growth of the initiator 
(a certain bond k), the avalanche lasts until the bonds 
selected to grow are those joining the growth interface 
after the growth of bond k (bonds "younger" than k). 
Note that all these bonds have the related random num- 
ber X smaller than the initiator one. If the bond selected 
by the dynamics was on the perimeter before the initia- 
tor growth, then the avalanche stops. In the asymptotic 
limit, due to the step shape of the histogram, only bonds 
with X < Xc grow. We call Q{s] x) the size distribution 
of avalanches whose initiator is associated with a num- 
ber equal to x. It is shown for bulk IP, both through nu- 
merical simulations and theoretical calculation, Q{s; x) is 
scale invariant (i.e. is a power law), only if the variable 
X of the initiator is equal to Xc- If x < Xc Q{s;x) has 
an exponential cut off at a typical size sq ~ {xc — x)~'^ 
with <T > 0. The avalanche distribution Q{s; Xc) for bulk 
IP without trapping has a power law shape with an ex- 
ponent r''"' - 1.60 ±0.02 |,|l|l. The dynamical activity 
near the boundary can be characterized by the distri- 
bution of the avalanches whose first bond (initiator) is 
located on the boundary. We performed a set of about 
10'^ numerical simulations of IP without trapping, of size 
3L X 5L with L = 128, lasting 4 x 10^ time steps and col- 
lected the statistics of boundary avalanches from the last 
2 X 10^ time steps, in order to ensure that the system 
is in its asymptotic critical state. To identify the sin- 
gle avalanche, we followed ||^, by adding the condition 
that the initiator of the avalanche is on the boundary. In 
FIG. we show the behaviour of the boundary avalanche 
distribution. We find: r'*"'" = 1.56 ± 0.05. This value is 
very near to the bulk value, and we can conclude from our 
numerical analysis that bulk and boundary avalanches 
have the same distribution. In section V we will derive 
analytically this result. 
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III. RUN TIME STATISTICS AND FIXED SCALE 
TRANSFORMATION 

In this section we introduce the theoretical tools we 
used to compute the boundary fractal dimension d'^"^ of 
the infinite IP cluster and the boundary avalanche expo- 
nent r""''. Our strategy combines Fixed Scale Transfor- 
mation (FST) and Run Time Statistics (RTS) gQ- 
We describe briefly the FST approach and we focus more 
on the RTS. 

FST is a lattice path integral scheme allowing one to 
evaluate the spatial correlation properties of the inter- 
section between an infinite fractal cluster and a straight 
line. This approach is based on the statistical invariance 
of the correlation properties under a parallel translation 
of the intersecting line (valid for fractals which are homo- 
geneous, at least in the translation direction). In particu- 
lar, it is possible to compute the probabilities Cq,Ci and 
C2 related to the configurations 0, 1,2 of the fine grain- 
ing process of FIG.|[ For the normalization condition it 
follows: 



Co + Ci + C2 = 1 



(3) 



From these probabilities one can compute the fractal di- 
mension of the intersection by: 



log{Co + Ci + 2C2) _ log{l + C2) 



log2 



log2 



(4) 



As usual, due to the intersection dimension rule, the 
fractal dimension D of the analysed cluster is given by 
D = 1 + d. The probabilities Co,Ci and C2 are com- 
puted through the statistical weights of growth paths, 
once a stochastic dynamical formulation of the model is 
given. This means that the use of FST is straightfor- 
ward whenever a simple calculation of the growth paths 
on the lattice is possible. In the present case, there are 
two problems to overcome in applying the FST. 

Firstly, the fractal properties of the system depend 
on the distance from the boundary (Cq = Co{z),Ci = 
Ci{z),C2 = C2{z)), this extrapolation from the inter- 
section dimension to the global dimension is no more al- 
lowed. Moreover, what we actually can compute with the 
FST method are the local (near to the boundary) cor- 
relations orthogonal to the boundary, while the fractal 
dimension of the intersection set parallel to the bound- 
ary is given by the correlation properties parallel to the 
boundary. However, since the crossover of the fractal 
dimension from the boundary to the bulk is very slow 
(logarithmic), one is allowed to assume that the cluster 
is "locally" isotropic. In this case transversal and hori- 
zontal correlations in a thin (with respect to the system 
size) strip parallel to the boundary share similar prop- 
erties. For the same reason, we can evaluate the fractal 
dimension d of the intersection between the cluster and 
a straight line parallel to the lateral boundary, through 
the first neighbors correlations orthogonal to the same 
boundary at the same distance. 



Secondly for IP (and for any other model with de- 
terministic extremal dynamics) the calculation of the 
growth paths is extremely difficult, because the weight 
of a path cannot be written as the product of the prob- 
abilities of the single steps composing it. The extremal 
dynamics of IP is deterministic, and the disorder appears 
only as a realization of quenched random variables. This 
implies that to evaluate the statistical weight of a given 
path we have to perform an average over all the quenched 
disorder and this average does not factorize itself in the 
product of the averages of the single steps composing the 
path. The latter problem is solved by the introduction of 
the RTS transformation. This transformation allows us 
to represent a quenched-extremal process like a stochas- 
tic dynamics. 

As regards the RTS (for a more detailed discussion see 
(P), the starting point is, at each time step t, to consider 
an effective probability density pi,t{x) for the random 
number Xi associated to each bond i of the growing in- 
terface dCt- This density depends on the growth history 
of the dynamics. In fact, pi^t{x)dx gives the probability 
that the variable Xi for the bond i at time t is in the 
interval [x,x + dx\, conditioned by the past growth dy- 
namics of the cluster. If a bond i does not belong to the 
cluster, or to the growth interface, its effective probabil- 
ity density is the flat one. Meanwhile, the bonds on the 
growth interface show a more interesting form of distri- 
bution. Once the densities Pi.t{x) for each bond i on the 
interface are known, one can calculate the growth prob- 
ability distribution {/i^,*} (i.e. the probability of being 
the minimum on the interface) at that time step for each 
interface bond (see appendix A): 



= / dxpi^t{x) TT / 



jeact-{i:} 



dypj,t{y) 



(5) 



where dCt — {i} represents the growth interface at time 
t except for the bond i. The effective probability den- 
sity of any surviving bond j at time t + 1 on the inter- 
face must then be updated, conditioned to the previous 
growth history at time t, i.e. the growth of the bond i. 
The corresponding equation is (see appendix A): 



Pi,t 



/ dypiAv) n / 



kedCt-{i,j} 



dzpk,tiz) 



(6) 



where dCt — {i, j} is the growth interface except for bonds 
i and j. New bonds added to the perimeter are assigned 
an effective probability density according to a uniform 
distribution in [0, 1], as no information is available about 
them till this time step. 

The above formalism allows us to write the statistical 
weight of a path as the product of the probabilities of 
individual steps. 
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IV. COMPUTATION OF THE BOUNDARY 
FRACTAL DIMENSION 

In order to combine the FST and the RTS approach, 
we need to have scale invariant growth rules (we want to 
compute a critical exponent, the fractal dimension, and 
the result cannot depend on the scale). The extremal 
dynamics of IP is known to be independent on the choice 
of the initial distribution of quenched variables. By using 
this symmetry, one can show that the scale invariant dy- 
namics for block variables is identical to the microscopic 
dynamics . The FST performs the computation of the 
correlation properties of a given structure by considering 
only the growth processes inside a growth column (FIG. 
^) . This approximation has been shown to be a good one 
for the dielectric breakdown model (DBM) jl^, and for 
bulk IP |]. Since eqs. (H), (^) involve all the variables on 
the perimeter of the growing cluster, a limitation of the 
process in the growth column destroys these correlations, 
leading to compact clusters H] . The solution to this prob- 
lem is given by observing that, as the critical avalanche 
size distribution is a power law, the statistical properties 
of a generic one (i.e. an avalanche whose initiator i has 
Xi — Xc) are then scale invariant. Then if one considers 
the dynamical evolution of a generic critical avalanche 
inside the growth column one obtains the scale invari- 
ant correlation properties (i.e. Cq, Ci and C2) needed 
to compute the fractal dimension. This can be done by 
modifying the eqs. (||), (^, in order to take account the 
dynamical evolution of a single critical avalanche. We 
consider a growth column on the perimeter of the infinite 
structure {t — >• 00). The starting point is the observation 
that scale invariant asymptotic avalanches begin with an 
initiator at x = Xc, due to the asymptotic shape of the 
histogram. All the memory of the past growth history is 
then contained in the requirement that the initiator has 
X — Xc- Then one is allowed to consider explicitly only 
the bonds grown after the growth of the initiator. 

The RTS dynamics corresponding to the local scale 
invariant dynamics, is obtained by 

• considering only bonds inside the growth column; 

• imposing that any "active" bond i in the column 
can grow only if the value of its variable Xi is less 
than Xc = 1/2. The idea is that if Xi > Xc for all 
the bonds in the growth column, growth will occur 
at some other place in the structure outside the 
growth column (it coincides with the definition of 
scale invariant avalanche); 

• imposing that the initial bond (i.e. the initiator), 
which is the largest of the variables participating 
to the growth process, has exactly Xi = Xc- 

In this way we modified the Eqs. (||) and (||), limiting the 
product over the perimeter variables to variables inside 
the growth column. 



Because of the presence of a lateral surface, this model 
is intrinsically anisotropic, and consequently we have to 
introduce some modification to the usual way of perform- 
ing the FST for the bulk IP. The anisotropy of the envi- 
ronment implies a breaking of symmetry in the FST basic 
configurations in FIG.^. Then, due to the presence of the 
boundary, the probabilities Cq and Ci are not equal in 
this case. 

Through the FST one may compute directly the matrix 
elements M,,- and from the relation: 




Moo Mio M20 
Mm Mil M21 
M02 M12 M22 




(7) 



it is possible to evaluate C2 and, by using Eq.d), d. In 
this case 



and 
C2 



Moi = Mio = 



M12MQ 



M12 + Af2l(Mo2 - M12) + Mi2{Mo2 - M22) 



(8) 



(9) 



The anisotropy of the environment is also introduced 
in the lateral boundary condition of the growth column 
where the FST calculation is performed. At the left side 
of the column we impose the presence of a rigid wall and 
at the right side the paths are allowed to go out and then 
to return inside the growth column, as can be seen in 



FIG.hO 



In this way we have obtained the results shown 



in Tab. Ill, where the fractal dimension for increasing or- 



der n (the path length) of the FST computation is given. 
We used a power law fit (FIG.^l]) to extrapolate df^^{n) 
to n = cx) and obtained ^""''(00) ~ 0.70. 

A similar approach has been applied to IP with site 
trapping, in particular: when a growth path produces a 
closed region surrounding the initial pair configuration 
Ci (see FIG. ^0|), it stops and its statistical weight con- 
tributes to the matrix elements M^^i, since the empty 
right (or left) site above the initial configuration cannot 
be occupied anymore. The results are shown in Tab. IV 
and in FIG.^ We have extrapolated our results to 
71 = CX3 by using the following function (see FIGjl^): 



d«;^'-(n) = -exp{-n°') 
n 



(10) 



with a — 0.66. This extrapolation gives df"''(oo) ~ 0.66. 



V. COMPUTATION OF THE BOUNDARY 
AVALANCHE EXPONENT 

We now propose a simple theoretical scheme for the an- 
alytical calculation of the boundary avalanche exponent 
^sur jp -^^jthout trapping, based on the RTS and the 
FST ideas, which has been successfully applied to bulk 
IP 1. 
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The following functional form for the avalanche size 
distribution is assumed: 



(11) 



J2Q(.s;x) = 1 Va; e [0,1], 



where xf^"^ = 1/2 is the critical threshold. The function 
/(x) has the following properties: lim^^^o fi^) — 
and for large values of x one has f{x) ~ e~^. Since the 
size s of the avalanche also includes the initiator, the 
normalization condition for the eq.(pT|) is: 



(12) 



Usually eq.(|ll|) holds for s ;» 1. However, if we consider 
the dynamics at a certain scale £, we can use eq.(^l|) to 
describe the statistics of avalanches at that scale. In the 
limit t oo, for x — x^"*", the asymptotic behaviour 
described by eq. (|ll]) holds for smaller and smaller values 
of s as ^ is increased. The deviations from the pure power 
law behaviour are integrated out into the dynamics at 
scale £. For £ ^ 1 we are allowed to suppose that eq. (|ll| ) 
holds from s = ltos>>l. In this case the normalized 
form of eq.(|0), for x = cc^"'' is: 



Q{s;xD^ 



(13) 



The denominator of eq. ( [T^ ) is the Riemann zeta fun 
tion, ((r''"'). 

From Eq.(|l3|), valid if the initiator is at Xc one has: 



Q(s = i;xr) 



1 



C(t''"'')' 



(14) 



To obtain an analytic estimation of the boundary 
avalanche exponent r'*"'' one has to evaluate the left hand 
side by taking into account the boundary conditions near 
the avalanche, together with the presence of the bound- 
ary. Then inverting eq. ( p^ it is possible to measure r^"'' . 
Let us evaluate Q{s = l;xf^^). The event s = 1 means 
that after the growth of the initiator with variable a;^"'' 
the avalanche stops. Thus, we consider the initiator i 
that grew at time tg E^^id we compute the probability 
that the avalanche stops time to + 1- This happens when 
all the descendant bonds of the initiator have variables 
larger than x^"''. In fact, if at least one descendant of i 
had variable lower than xf^^, the avalanche would con- 
tinue because this variable would be the minimum one 
on the whole perimeter. In order to evaluate this prob- 
ability we need to take into account the environment of 
the initiator. In FIGjlJ a-b we schematize all the possi- 
ble boundary conditions for the initiator bond. We con- 
sider only the nearest neighbours of the initiator because, 
asymptotically, the avalanches on the perimeter are influ- 
enced only by the environment near the zone where the 
avalanche evolves. That is, they are affected by other 
branches of the aggregate which have some perimeter 



bonds affected by the avalanche. The presence of the 
boundary is implemented by allowing only the right and 
the vertical bond to grow in FIG. ^ 

For these two cases we can evaluate the probability 
that the avalanche stops immediately after the growth 
of the initiator, conditioned by the assigned boundary 
conditions. The exact value of this probability is given 
by the average of the two cases. In order to calculate 
the statistical weights of configurations (a), (b) and (c) 
of FIGjl^ we use the void distribution P(A) of the ran- 
dom anisotropic Cantor set whose generators have local 
(near the boundary) probabilities z = 0,1,2 given 
by the FST calculations performed in the previous sec- 
tion. We are then allowed to use P{X) with the weights 
obtained by FST because for IP the perimeter has the 
same statistical properties as the bulk of the structure. 
Obviously, the void distribution we obtain is a local one, 
since the probabilities Ci for the Cantor set orthogonal 
to the boundary are dependent on the distance from it. 
In practice, only the P(A = 0) can be computed with a 
reasonable degree of accuracy, because it depends only 
upon the local properties of the set. When the size A 
of the void is not small with respect to the system size, 
the implicit assumption that the Ci are independent of z 
becomes inconsistent. 

We report the expression of P(A = 0) from ||l^ in 
terms of C2 and Ci : 

The weight of configuration (a) is: 

yy(a) = 1 _ p(A = 0). 

The weight of configuration (b) is: 

W'-''^ = P{X = 0), 

The fixed point values of C2 and Ci obtained from FST 
calculation ofd'pgrp in the previous section are C2 — 0.628 
and Ci ~ 0.249. If we introduce these values in eq.([l5|) 
we get: P(A = 0) ~ 0.648. The probability Q{s = l;Xc) 
to have an avalanche of duration s = 1 is 

Q{s = l;xc) = (1 - x,f{l - P(A = 0))+ 

(1 -Xc)P(A = 0) = 0.412 (18) 

At this point, in order to find r*"'" we should solve the 
equation: 



(16) 



(17) 



0.412 



1 



The numerical solution of eq.(19) gives 



1.55 



(19) 



(20) 



in very good agreement with our numerical findings. The 
above scheme is, however, too simplified to account for 
trapping. In fact, the method is based on the first growth 
step inside an avalanche, while trapping becomes relevant 
at higher orders (see Table III and Table IV). 
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VI. CONCLUSIONS 

In this paper, we have presented, in analogy with usual 
critical phenomena, the study of boundary effects in in- 
vasion percolation with and without trapping. Near a 
boundary one deals with a qualitatively different rate 
of occupation. This is reflected in a lower fractal di- 
mension of this part of the cluster. Numerical simula- 
tions give surface fractal dimensions d^^^ = 0.65 ± 0.02 
and df;^'' = 0.59 ± 0.03 for IP without trapping and IP 
with site trapping respectively. These two values are 
smaller than the bulk values. Meanwhile, simulations 
for the asymptotic shape of the histogram distribution 
and for the boundary avalanche distribution for IP with- 
out trapping, show that the boundary does not affect 
these quantities. The histogram self-organises into a 
theta function with threshold x^"'' = 1/2 and the bound- 
ary avalanche distribution is characterized by an expo- 



Prob{A\B) 



Prob{A fl B) 
Proh{B) 



(A3) 



The events A and B are respectively A = {x < Xj < 
X + dx) and B = {xi = min^ggc^ Xk)- By definition of 
"effective probability density" , we can write: 

Prob{t + l;x < Xj < a; -f dx) — dx pj,f+i{x). (A4) 

However, using conditional probability, we can also write 

Prob{t -\- I; X < Xj < x + dx) = 



Prob{t; {x <Xj < x + dx) \ {xi — min Xk)) ~ 

k£dCt 

Probjt; {x < Xj < X + dx)f]{xi = minfcggct Xk)) 
Prob{t; Xi = minkedCt ^k) 



(A5) 



nent r™"- = 1.56 ± 0.05, very near to the bulk value The numerator of ^ can be written as: 



^buik _ i,go±0.02. We are also able to present a theoret- 
ical scheme to compute analytically the relevant critical 
exponents d*"'' (for both IP with and without trapping) 
and r''"'' (for IP without trapping only) near the bound- 
ary. Our theoretical results d*""" ~ 0.70, d^^^ ~ 0.66 and 
^sur ^ 2^ good agreement with the numerical 

data. Authors acknowledge S. Cornell for suggestions, 
GC acknowledge the support of EPSRC. 



APPENDIX A: DERIVATION OF THE RTS 
EQUATIONS 

In Invasion Percolation a bond grows at time t if its 
variable is the minimum one at that time. Then we can 
write: 



Prob(t; {x < Xi < x + dx) [^(^ 



min Xk)) 
kedCt 



Prob(t; {x < Xj < x + dx) ^{xi — mm Xk)) 



^ dx pj^t{x) dy pi^tiy) Y\. [ dupk,t{u) 

(A6) 

This gives the probability that, at time t, x < Xj < 
X + dx, and at the same time xt = min^ggct ^k) (i.e. 
Xi = y Cz [0, x] and for all the other k G dCt, Xk > y. The 
denominator of the right term in Eq.(A6) is simply pi^f 
Then we have: 



Pj.t+i{x) 



pjAx) r 



/ dyp,^{y) \Y / dzpk,t{z) 



dx Pi,tix) 

kedCt-ii} 



PkAv)dy 



(Al) 



This gives the probability that, at time x < xi < x+dx 
and at the same time Xi is the minimum on dCt (i.e. that 
every other bond variable in dCt is between x and 1). 
By integrating eq.(Al) one can finally write the growth 
probability pi,t for the bond i at time t [|l^,p|: 



Pi,t = Prob{t;Xi 



min Xk) 

kedCt 



/ dxpi^x) TT / 



Pk.t{y)dy 



(A2) 



keact-li} ' 



To update the effective densities pj,t{x) of generic bond 
not grown j £ dCt to obtain pj,t+i{x), we make use of 
the law of conditional probability: 
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L = 64 


L = 128 


L = 192 


L = 256 


H = 2 


0.66 ± 0.02 


0.64 ± 0.02 


0.63 ± 0.02 


0.62 ± 0.02 


H = 3 


0.70 ± 0.02 


0.66 ± 0.02 


0.64 ± 0.02 


0.63 ± 0.02 


H = 4 


0.73 ± 0.02 


0.68 ± 0.02 


0.66 ± 0.02 


0.65 ± 0.02 



TABLE I. Boundary fractal dimension of IP without trap- 
ping for different system sizes L and different values of H. 
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0.59 ±0.03 


0.58 ±0.03 


H = 3 
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0.62 ± 0.03 


0.60 ± 0.03 


H = 4 
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TABLE II. Boundary fractal dimension of IP with trapping 
for different system sizes L and different values of H. 
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n 


2 


3 


4 


5 


6 


7 




00 




0.453 


0.576 


0.622 


0.641 


0.651 


0.657 




0.664 



TABLE IV. Values of the boundary fractal dimension of IP 
with site trapping with respect to the order n of computation. 



n 


2 


3 


4 


5 


6 


7 








O.l."):! 


{).57() 


().():!2 


{).()57 


O.GTi 


().()79 




0.702 



TABLE III. Values of the boundary fractal dimension with 
respect to the order n of computation, for IP without trap- 
ping. 
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invader 

FIG. 1. Set up of numerical simulations. An invading (not 
yet percolating), cluster is shown. Only the bottom left part 
of the cluster will be considered for the statistics. 



-0.3 
-0.5 — 



FIG. 4. Collapse plot of [d(«) -d'""']Zog(L) for the different 

sizes, in log — linear scale, for IP without trapping {H = 4). 
z and L are in lattice units, then the normalized distance is 
dimensionless. 



0.90 




FIG. 2. This picture shows the entire cluster. The region z/L from the boundary {H = 4). z and L are in lattice units, 
of interest in which statistics is taken is the lower-left one. then the normahzed distance is dimensionless. 
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L=64 




L=128 
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L=192 




L=256 


1 

1 





0.0 0.2 0.4 0.6 0.8 1.0 

z/L 



FIG. 3. Behavior of the fractal dimension of the intersec- 
tion sot of IP without trapping versus the normalized distance 
z/L from the boundary {H = 4). z and L are in lattice units, 
then the normalized distance is dimensionless. 




0.0 0.2 0.4 0.6 0.8 1.0 



FIG. 6. The histogram distribution of bulk perimeter vari- 
ables in IP without trapping is compared with the distribution 
of variables near the boundary of the system, after 5 x 10^ 
time steps. The two distributions coincide. 
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» border av. distr. 
t(s)=a»s^, 1=1^6, A'C=0.05 
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FIG. 7. Border avalanche distribution of IP without trap- 
ping in log — log scale. The least square fit gives a slope 
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FIG. 10. A possible path of growth at the G*** order. The 
invasion proceeds along the arrows, from one black point to 
another one. The number near each arrow is the growth time. 





FIG. 9. The growth column used in the FST scheme. The 
thick bond is the initiator of the avalanche, with x = Xc, and 
the dashed bonds are the bonds of the perimeter after the 
initiator's growth. 
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1/n expC-n"**) 

FIG. 12. Extrapolation of the FST fractal dimension 
d|"''for IP with site trapping. 
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(a) (b) 

FIG. 13. (a)-(b): Boundary conditions for a boundary 
avalanche after the growth of the initiator. ( (•) indicates the 
cluster sites and (o) the perimeter ones: the filled segments 
are grown bonds and the dotted ones are the descendant of 
the initiator. The left boundary is shown. Its effect is to re- 
duce the maximum number of perimeter bonds in which the 
avalanche can go (2), with respect to the bulk (3). 
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